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Abstract 

We report a Molecular Dynamics study of homogenous bubble nucleation in a Lennard-Jones 
fluid. The rate of bubble nucleation is estimated using forward-flux sampling (FFS). We flnd that 
cavitation starts with compact bubbles rather than with ramified structures as had been suggested 
by Shen and Debenedetti ( J. Chem. Phys. Ill, 3581 (1999)). Our estimate of the bubble- 
nucleation rate is higher than predicted on the basis of Classical Nucleation Theory (CNT). Our 
simulations show that local temperature fluctuations correlate strongly with subsequent bubble 
formation - this mechanism is not taken into account in CNT. 
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1. INTRODUCTION 



Anyone who has ever sprinkled water droplets into a pan with hot frying oil is familiar 
with the phenomenon of explosive boiling: the liquid water will be heated far above its 
equilibrium boiling point before the vapor phase nucleates suddenly and violently. As this 
example illustrates, explosive boiling is very common, and usually undesirable in practical 
situations. Bubble nucleation is relevant in many different contexts: in addition to explosive 
boiling^, it plays a role in phenomena as diverse as cavitation erosion^ and sonochemistry^, 
and it affects the design of high-efficiency heat exchangers^. In spite of the practical relevance 
of these phenomena, the mechanism by which the vapor phase nucleates from a homogeneous, 
super-heated liquid is still under debate. 

The standard theory used to describe the bubble nucleation phenomenon is the so-called 
Classical Nucleation Theory (CNT). CNT is commonly used to predict the rate of nucleation 
and estimate the height of the free-energy barrier.-'^'^'^i^ CNT assumes that the growing 
droplet of the stable phase within the metastable one is characterized by its bulk properties, 
it is incompressible, and it has a spherical shape. Moreover, the liquid-vapor interfacial 
free-energy does not depend on the curvature of the droplet, condition that is satisfied at 
the surface tension, that coincides with the equimolar dividing surface in case of absence of 
Gibbs absorption. 

From an experimental point of view, CNT is in general used to predict the nucleation 
rate, both in bubble nucleation and in vapor condensation experiments. Early measurements 
of vapor condensation performed by Blander and Kat^iS turned out to be consistent with 
the CNT predictions^'ii'i^. However, more recent measurements using a variety of different 
techniques and niaterials^'^'^>^'^>^'^'^'^'2^'^'2ii2^>2^i22i22ii22 demonstrated that although 
CNT correctly predicts the isothermal supersaturation-dependence of the bubble-nucleation 
rate, it fails to predict the correct temperature dependence: CNT predicts nucleation rates 
that are too low at the low temperatures and too large at high temperatures (even allowing 
for the fact that error bars may span several orders of magnitude). To our knowledge, 
there exist no experimental techniques that can directly probe the rate-limiting step in the 
nucleation process, namely to formation of "critical" bubbles. 

In an attempt to improve upon the predictions of CNT, Zeng and Oxtoby^ have used 
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Density-Functional Theory (DFT) to arrive at an estimate of bubble nucleation rate. In 
particular, Zeng and Oxtoby studied the liquid-to-vapor phase transition in a Lennard- 
Jones fluid and found that CNT underestimated the nucleation rate by more than 15 orders 
of magnitude. In order to explain the discrepancy between the DFT predictions and the 
CNT estimates, Oxtoby and co-workers^ii^ pointed out that CNT neglects all curvature 
corrections to the surface free-energy of a droplet. Another curious feature of CNT that was 
already pointed out by Cahn and Hilliard^'^ and discussed in some detail by Oxtoby and 
Evans^^ is the fact that this theory predicts a finite nucleation barrier at the point where 
the metastable liquid undergoes spinodal decomposition — this topic is still the subject 
of much debate.^i22 Delale et al.~^ proposed a phenomenological estimate of the minimum 



wor. 
ref. 



i of bubble formation, using experimental data on super-heating. Using this approach. 



38l estimated the steady-state bubble nucleation rate and found that, as with the CNT 



estimates, there was a large discrepancy between the predicted computed and measured 
bubble-nucleation rates. 

In 1999, Shen and Debenedetti^ reported a Monte Carlo (MC) study to estimate the 
free-energy barrier for bubble nucleation in a monatomic fluid of particles interacting via a 
truncated Lennard- Jones potential. These authors computed the free-energy barrier using 
Umbrella Sampling^"^, and monitored the progress of bubble formation using the total 
density of the system as an order parameter. A geometric analysis of the resulting void 
structures in the superheated liquid revealed that the formation of the vapor phase pro- 
ceeded via the formation of ramifled, system-spanning void structures. This observation 
disagrees qualitatively with the predictions of CNT and DFT, both of which assume that 
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the critical nucleus is a compact, spherical bubble. It should however be noted that ref. 
studied a system very close to the point where the system is expected to undergo spinodal 
decomposition. 

The aim of the present paper is to study the kinetics of bubble nucleation under conditions 
similar to those studied in ref. l39|. In order to study the time evolution of bubble nucleation, 
we must use Molecular Dynamics (MD) simulations and because it is then better to use 
an interaction potential such that the forces vanish continuously at the cut off, we could 
not use exactly the same model as Shen and Debenedetti. Rather, we studied a system 
of particles interacting via a truncated, force-shifted Lennard- Jones potential. However, 
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as we shall argue below, we study thermodynamic conditions that are comparable. Using 
Forward-Flux Sampling (FFS)^"^, we compute the bubble nucleation rate and compare our 
numerical estimate with the prediction of CNT. FFS is well suited to simulate nucleation 
processes because it does not assume that the nucleation process is slow compared to the 
time it takes all other degrees of freedom to equilibrate. The latter assumption is implicit 
in the Umbrella Sampling scheme. 

The rest of this paper is organized as follows: we first describe the simulations techniques 
used in our work, then we present our results and discuss the imphcations. 

2. SIMULATION METHODS 

In order to probe the mechanism of bubble nucleation, we simulated a system consisting 
of 3375 particles interacting via a truncated and force-shifted (TSF) Lennard- Jones (LJ) 
pair potential: 

(r - Tc) . (1) 

The TSF-L J potential is a modification of the 12-6 Lennard- Jones potential. The second and 
third terms on the right hand side of Eq. [T]have been chosen such that both the potential and 
its first derivative vanish continuously at r^. = 2.5a (where a is the particle diameter). The 
reason why we use the TSF-LJ rather than the full LJ potential is that during nucleation 
the system becomes inhomogeneous. Under those conditions it is much simpler to use a 
finite-ranged potential than an interaction with an r^^-tail. 

In what follows, we use reduced units: all physical variables are expressed in terms of 
the Lennard- Jones units a, m, e, where a is a measure for the diameter of a Lennard- Jones 
particle, m denotes its mass and e is the depth of the attractive well of the (untruncated) 
LJ potential. 

We start our simulations by equilibrating the system in an NPT ensemble by means 
of MD simulations with a Nose-Hoover thermostat^i^i^ and an Andersen barostat^'' (Ap- 
pendix A). In order to study the bubble nucleation phenomenon, we use MD rather than 
MC simulations, because standard MC algorithms assume that the system is always locally 
in thermal equilibrium - an assumption that may not be justified in the present case. In all 
simulations, we adopt cubic periodic boundary conditions. The equations of motion are in- 



VTsrir) = VLj{r) - VLjir^) - ^^^^ 

ar 



5 



tegrated using a leap-frog algorithm^ with an integration time step of At* = 0.00046, which 
corresponds to a time step of about 1 femto-second (using the Lennard- Jones parameters 
for Argon). In order to vahdate the MD code used in the present work, we verified that it 
could successfully reproduce the liquid-vapor coexistence curve of the TSF-LJ model that 
was computed by Errington et al.— (Appendix B). 

In an NPT ensemble the simulation box adjusts itself in order to keep both the temper- 
ature and the pressure of the system constant. When a system undergoes a hquid-to- vapor 
during an NPT simulation, the simulation box tends to expand dramatically: in our simu- 
lation runs we limit such "explosions" by interrupting the runs before the liquid phase has 
evaporated completely. 

In order to study bubble nucleation, the system should be large enough to accommodate 
a "critical" bubble (i.e. a bubble that is equally likely to grow or shrink). On the basis of 
CNT we expect that the size of the critical bubble is inversely proportional to the degree of 
super-saturation. If we carry out simulations at a low supersaturation, the critical bubble 
size can become very large and this would necessitate the use of very large simulation boxes 
(containing a large number of particles). On the other hand, if we choose the supersaturation 
too large, the liquid- vapor phase transformation will not proceed via nucleation and growth 
but through spinodal decomposition. These two constraints put rather narrow limits on the 
values of supersaturation that can conveniently be studied. 

In order to obtain a rough estimate of the critical bubble size, we make the CNT assump- 
tion that the chemical potential of the vapor inside the critical bubble is the same as that 
of the liquid outside. If we use the fact that the molar volume of the vapor phase is much 
larger than that of the liquid phase, this implies that pressure inside the critical bubble is 
the same as the equilibrium vapor pressure at the imposed temperature. The radius of the 
critical bubble can then be approximated as: 

"■^"^ " P,p{T)-P ' 

where 7 is the liquid-vapor surface tension, P is the imposed external pressure of the liquid 
and Pyp{T) is the equilibrium vapor pressure at the temperature T. 

Shen and Debenedetti^^ studied bubble formation in the truncated Lennard- Jones system 
at a saturated vapor pressure of P* = 0.046 and various degrees of super- heating, defined 
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as S* = (T* — Tsat)/Tsati where Tgat the vapor-hquid coexistence temperature at the given 
pressure). For 5=9%, CNT predicts a free-energy barrier of 37/c bT for bubble nucleation.— 



In the present study, we do not use exactly the same model as ref. 



39, but we do use the same 



degree of super-heating (5* = 9%). Although the FTS-LJ model is not identical to the model 



used in ref. 
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I, our simulations were carried out at approximately the same distance from 
the critical point (T*/T^*=0.914). For this reduced temperature and supersaturation, we 
expect that a system of 3375 particles can accommodate a critical bubble. To summarize: we 
carried out simulations at a temperature T* = 0.855 and a reduced pressure P* = 0.026. The 
coexistence temperature at this reduced pressure is T*^^ = 0.785. The critical temperature 
of the TSF-LJ model was estimated by Errington et al.^- to be T* = 0.935. At coexistence 
i'^sat = 0.785, P*^^ = 0.026), the number densities of the coexisting liquid and vapor are, 
respectively, p^coex = 0.668 and plcoex = 0.043. The number density of the superheated 
hquid at T* = 0.855 and P* = 0.026 is = 0.580. 



2.1. Liquid- vapor surface tension 

In order to compare our numerical estimate for the nucleation rate with the prediction 
based on CNT, we need to know the value of the surface tension 7. As both the critical 
temperature and the surface tension are sensitive to the long-range part of the intermolecular 
potential, it is important to have a reliable estimate of 7 for the TSF-LJ potential truncated 
at Tc = 2.5(7 for the state point studied in the simulations. 

To estimate 7, we have performed MD simulations of the TSF-LJ liquid- vapor coexistent 
systems at T* = 0.785 (P* = 0.026) and at T* = 0.855 {P* = 0.046) separately. We 
obtain the surface tensions of the planar interfaces: 7* = 0.204 ± 0.007 at T* = 0.785, and 
7* = 0.098 ± 0.008 at T* = 0.855 (Appendix C). Lutsko has used DFT to estimate the 
liquid-vapor surface tension of the TSF-LJ system at the same state points and obtained 
Yj^p^ = 0.198 at T* = 0.785, and Ydft = 0-119 at T* = 0.855.-^'^ These values for the 



39. In view of 



surface tension are very similar to the ones for the model studied in ref. 
the statistical errors in the simulation data, the agreement between the DFT estimates for 
7* and the simulation results is quite good. We did not compute the surface tension of a 
super-heated system at T* = 0.855 and P* = 0.026 because, under this thermodynamic 
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condition a flat interface is not tliermodynamically stable. 



2.2. The order parameter 

Even in the stable liquid phase, spontaneous density fluctuations occur. In a superheated 
liquid, some of these local density fluctuations may grow to nucleate a bubble of the stable 
vapor phase. In order to follow the nucleation process we need an "order parameter" that 
provides a quantitative measure of the degree of transformation from the liquid to the vapor 
phase. In the present simulations we used the volume of the biggest bubble as an order 
parameter. This is a local order parameter that is sensitive to the formation of large, 
connected void spaces. The present deflnition of the order parameter is different from the 
one used by Shen and Debenedetti, who used the overall density of the system as a global 
order parameter. We will return later to this difference in the choice of order parameters. 

We used the following procedure to identify the largest connected low-density region in 
the superheated liquid: 



1. Following the method used in ref. l52|, we construct a three-dimensional grid with cell 
mesh of (O.Scx)^ inside the simulation box [a is the LJ diameter). 



2. As in ref. |53|, we use Stillinger's cluster criterion^^ to identify the nearest neighbors of 
each particle: according to our deflnition, two particles are neighbors if their distance 
is less than a cut-off = 1.6cr, corresponding to the flrst minimum of the radial 
distribution function of the liquid at the same thermodynamic conditions. Using this 
criterion, we compute the probability distribution of each particle's neighbors in both 
the liquid and the vapor phases. 

3. We then identify a particle as liquid-like if it has more than flve neighbors. We used 
this particular threshold, as essentially all particles in the homogeneous vapor phase 
have fewer neighbors and almost all particles in a homogeneous liquid phase have more 
(see Figured]). 

4. We deflne a spherical volume with radius rc around every liquid-like particle: all grid 
cells which are inside this spherical volume are labelled as liquid-like. The remaining 
grid cells are labelled vapor-like. 
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5. Finally, we perform a cluster-analysis on the vapor-like grid cells: a set of connected, 
vapor-like cells defines a bubble. We then select the largest of those and define this as 
the largest bubble. The local order parameter is the volume Wb of the largest bubble 
expressed in units of (for practical reasons, we multiply the bubble volume by 10 
and express it as an integer in FFS computations). 

2.3. Forward Flux Sampling with Molecular Dynamics 

FFS is a technique that allows us to track the pathways and flux of a rare phase tran- 
sition between two regions of the phase space, the liquid (Lug) and the vapor(Vuap). As 
bubble nucleation is a rare-event phenomenon, we incorporated FFS technique^"^ in the 
MD simulations of the super-heated fluid. There are specific complications that arise when 
FFS is combined with MD and, to our knowledge, this had not been attempted before. 

Once we have defined a parameter A that measures the progress of the nucleation process, 
we can define a value of A, denoted as Aq, that define the boundary of the metastable liquid 
region. By this we mean a value of A that is large enough to ensure that the overwhelming 
majority of the states of the meta-stable liquid have a lower value of A, yet small enough 
to ensure that occasional a spontaneous fluctuation will cross Aq. We stress that the choice 
of Aq is to some extent arbitrary and that our results do not depend on the precise choice. 
Similarly, we define a maximal value of A, denoted as A„, such that almost all trajectories 
that reach this value of A will continue in the direction of progressive evaporation. We 
denote the domain of the meta-stable liquid (A < Aq) by Lug and that of the vapor (A > A^) 
by Kap- 

As the transition from Lug to V^ap is rare, the system will initially spend most of its time 
in the L^g-basin with A < Aq. The key idea of FFS is to use a series of non-intersecting hyper- 
surfaces (interfaces) in phase space, each one identified by a value of the order parameter 
A > Ao, to drive the system from Lug to V^ap in a ratchet-like manner. 

The first step in the simulation is to compute the flux of trajectories that cross Aq in 
the direction of increasing bubble size, per unit time and per unit volume. We store the 
configurations of the system at the time that a "forward" trajectory crosses Aq. Once a set 
of such configurations have been collected at Aq, we fire off stochastic trajectories that start 
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from Aq and either reach the next interface Ai or fall all the way back to the initial basin 
Liiq (Ao). We then compute P(Ao|Ai), the probability that a trajectory that had crossed Aq 
continues to Ai, rather than return to the initial basin Lug. We collect the configurations 
at Ai where trajectories from Aq first cross. Next we use these configurations to start new 
stochastic trajectories that will either reach A2 or fall back to Aq. In this way we compute 
P(Ai|A2), the probability a trajectory started at Ai reaches A2 before returning to the initial 
basin Lug. We then iterate this procedure all the way to the final boundary A„. The FFS 
estimate for the nucleation rate is then expressed as the product of two terms: 

n-l 

RL,^,v.., = <^-UPim+i) ■ (3) 

1=0 

The above description of the FFS technique makes it clear that there is a problem when 
combining FFS with MD: FFS that the underlying dynamics is stochastic, such that differ- 
ent trajectories can be generated from the same initial point in phase space. In contrast, 
standard Newtonian dynamics is deterministic, hence there is only one trajectory emanat- 
ing from a given point in phase space. In order to combine FFS with MD simulations 
in an NPT ensemble, we make use of a weakly stochastic variant of the (deterministic) 
Nose-Hoover thermostat. To this end, we add a Lowe- Andersen (LA) thermostat to the 
Nose-Hoover thermostat^. The LA thermostat is a momentum conserving and Galilean 
invariant analog of the Andersen thermostat: it re-scales the relative velocity between two 
particles from a Gaussian distribution with a probability At*!/*, where At* is the MD time 
step and u* a tunable parameter indicating the strength of the coupling between the system 
and the thermostat (the smaller the z/*, the weaker the coupling). Therefore the dynamical 
properties of the system can be tuned by varying u*. The larger u* the more stochastic the 
dynamics of our system. Of course, we wish to ensure that the stochastic nature of the LA 
thermostat has only a minor effect on the djTiamical properties of the system. As a test, we 
computed how the value of u* affectes the self-diffusion coefficient in the FST-LJ system (see 
Figure [2]). We find that when u* < 100, the diffusion of the system with the LA thermostat 
is indistinguishable from the one obtained with the deterministic Nose-Hoover thermostat. 
In fact, there is little change in the diffusion coefficient as long as u* < 1000. However, when 
u* > 10000, the LA thermostat clearly affects the dynamics. 

In order to balance the added stochasticity to the need of keeping the system's dynamics 
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as Newtonian as possible, we choose v* = 500. As Figure E] shows, the added stochasticity 
is enough to ensure divergence of the trajectories in phase space during the time it takes to 
travel from one interface to the next.— 



3. RESULTS AND DISCUSSION 

To perform an MD-FFS of bubble nucleation in a superheated liquid system at T* = 0.855 
and P* = 0.026, we locate the boundary of the metastable liquid basin at Wh = 25, and 
set the interfaces at monotonically increasing values of the order parameter Wb- Wh^i = 37, 
Wb,2 = 57, Wb,3 = 92, Wb,4 = 120, Wb,5 = 180, Wb,e = 240, W^j = 300, Wb,s = 400, 
Wb,9 = 550, Wb,io = 700, Wb,ii = 820, 1^6,12 = 950, Wb,i3 = 1200, Wb,u = 1900. At each 
interface we store O(IO^) configurations, select 30-50 of these at random and fire O(IO^) 
trajectories from each of these. At a value of Wb ~ 950, approximately 50% of all trajectories 
continue to the vapor phase. This value of Wb then constitutes our operational definition of 
the critical bubble size. However, we stress that, as Wb is not a perfect reaction coordinate, 
the committor of individual configurations with Wb = 950 may differ from 50%. We have 
not attempted to optimize the choice of the reaction coordinate to make it coincide more 
closely with an equi-committor surface. At Wb =3500, the committor to reach the vapor 
phase is very nearly equal to one: we use this value of Wb as the boundary of the vapor 
phase. Once simulations have reached this point, they are terminated. 

Using FFS, we obtain the following estimate for the bubble nucleation rate at T* = 0.855 
and P* = 0.026: R = 9 x lO^^^^^cr^'^r^^. Using the Lennard- Jones parameters of Argon, 
this rate becomes R = lO^^cm^^s^^. The computed bubble-nucleation rate is much higher 
than the one measured in early experiments for a number of volatile organic liquids: as 
an example, the bubble nucleation rate from ref. [lO is 10^ — 10^ bubbles cm~^s~^ at a 



temperature about 0.89Tc and an ambient pressure 1 atm. 
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3.1. Comparison with Classical Nucleation Theory 



A CNT-based theoretical estimate of the bubble nucleation rate was given by Katz (Eq. 12 



in Ref . 



id): 



■CNT 



N 



27 



.1/2 



nniB 



exp{-pAG) 



(4) 



where the prefactor (A^ 



27 



Zeldovitch factor 



27 



1/2 



1/2, 



contains the number density of the liquid (A^) and the 
m is the mass of a molecule, 7 the liquid- vapor surface tension, 
and B a term that takes into account the mechanical equilibrium of the bubble: in cavitation 
experiments, B is always equal to one. /3AG is the free-energy barrier to form the critical 
bubble, equal to g^p^^^p^p for a spherical bubble. 

The liquid-vapor surface tension of a planar interface of the TSF-LJ fluid at T* = 0.855 
is 7* = 0.098 ± 0.008 from MD simulations. This leads to a free-energy barrier of (3 AG = 
47 ± 11. Computing the bubble nucleation rate by means of Eq. HI we find that CNT 
predicts R ~ 10~^^o"^'^r^^, with an estimated error of five orders of magnitude in either 
direction. The CNT prediction is six orders of magnitude less than the simulation results. 
However, in view of the statistical errors in the data, this discrepancy is only marginal. 
Using the DFT estimate for the surface tension 7* = 0.119 at T* = 0.855, we obtain a CNT 
estimate R ~ O(10^^^)(T^'^r^^, which is 21 ± 1 orders of magnitude smaller than the MD- 
FFS computation. Our finding is consistent with earlier theoretical studies that found that 
CNT estimates of bubble-nucleation rates were always lower than the DFT results^'^'^. 
Of course, the DFT estimate of the surface tension at T* = 0.855 refers to the value for a 
planar interface at coexistence. 



3.2. Bubble nucleation pathways 

Figure H] shows a typical example of a computed pathway for bubble nucleation in a 
super-heated liquid at T* = 0.855 and P* = 0.026. As can be seen from Figure HJ the 
growing bubble is compact, even the critical bubble. In this respect, our findings appear 
to differ qualitatively from the findings of ref. Ssl. As the present simulations were carried 
out at virtually the same distance from the critical point and at the same superheating, it 
is likely that the difference in the results is due to a difference in the techniques that were 
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used to follow the progress of bubble nucleation. In particular, in the present study we 
made use of a local order parameter (viz. the size of the largest compact bubble), whereas 



ref. 



39l computed the free-energy as a function of a global order parameter (viz. the density 



of the entire system). These two methods need not probe the same process. As discussed 



in ref. 



571 there are two ways in which a system can accommodate an increase in a global 



nucleation order parameter. The first is via the growth of a compact nucleus, i.e. the process 
that we normally associate with nucleation. The second is via the generation of many small 
nuclei that, together, yield the same value of the global order parameter. If the system 
is large enough and the surface tension small enough, the second route will be favored for 
entropic reasons. As the global order parameter increases, the many small nuclei may form a 



39 . The fact that in our FFS simulation 



percolating structure, as was indeed observed in ref. 
we do not observe this scenario but the nucleation of a compact bubble indicates that the 
latter pathway is kinetically favored. 

We stress that the use of an order parameter that is based on the size of the largest bubble 
would still allow us to observe the formation of ramified clusters if this were the favored route: 
it would simply show up as sudden changes in the size of the largest bubble as it merges 
with other (smaller) bubbles. Interestingly, we occasionally do observe such processes, be 
it on a small scale: we find that sometimes two small bubbles may merge to form a larger 
bubble - we also observe the converse process where a larger bubble spontaneously breaks up 
into two smaller bubbles. When this happens the volume Wb of the largest bubble changes 
discontinuously (see Figure [5]). However, unlike ref. Q we find that the resulting bubbles 
are always compact. 

The conclusion must therefore be that FFS would in principle allow for the formation 



of ramified structures, as observed in ref. 



39l. The fact that we do not see such structures 



suggests that this pathway is kinetically unfavorable. Such kinetic effects could not be 



39. There are, however, other 



observed in the equilibrium umbrella-sampling study of ref. 
examples of nucleation events where the system appears to prefer a higher free-energy route 
for purely kinetic reasons.— 



13 



3.3. Local Kinetic Energy of Bubble Nucleus 



There exists and extensive literature on the study of rare events in condensed phases 



(see, for instance refs. 



59 



60l ). In the numerical study of rare events, such as activated steps 



in chemical reactions or nucleation events, the reaction coordinate is usually constructed in 
configuration space. In other words: the most common situation is one where the success of 
a barrier-crossing trajectory does not depend on the initial momenta of the particles involved 
in the rare event. However, in bubble nucleation the situation may be different. To see this, 
consider the converse of bubble nucleation, viz. bubble collapse. Bubble collapse tends to 
proceed sufficiently rapidly that heat cannot be transported away quickly enough to keep 
the process at constant temperature. Hence, at the end of a bubble collapse the system may 
be locally very hot. This is the physical origin of sono-luminescence (see, e.g. ref. l6ll ). It is 
therefore reasonable to assume that the reverse process, bubble nucleation, will preferentially 
originate in a part of the system that is locally hotter than its surroundings. Such an effect 
cannot be studied using a technique such as umbrella sampling, as this approach is based on 
the assumption of local thermal equilibrium. However, FFS makes no such assumptions and 
therefore allows us to see if local temperature fluctuations affect the propensity for bubble 
nucleation. 

To investigate the relation between local heating and subsequent nucleation, we studied 
the FFS trajectories of a large number of successful bubble-nucleation events. We stress that 
the initial conditions in these simulations were not biased: i.e. we did not modify the initial 
temperature. Once a bubble had successfully nucleated, we traced back the trajectory to the 
metastable liquid basin and analyzed the local temperature along this trajectory. To improve 
the statistics, we did this for a number of other successful trajectories that started from the 
same initial conditions. There are several ways to define the relevant local temperature. We 
used the following procedure: we identified all particles in the vapor phase inside the largest 
bubble at each interface along FFS trajectory and then tracked the kinetic energy of those 
particles from the basin to the end of FFS trajectories during the bubble nucleation process. 
As an example, we show the ensemble of the configurations collected at each interface from 
the successful pathways at T* = 0.855 and P* = 0.026 (T^^^ = 0.785) in Tabled 

In order to see if temperature has an effect on the propensity for bubble nucleation we 
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compared the local temperatures of bubbles that were committed to grow with that of all 
bubbles that have reached the same size. We can then compare these two temperatures with 
the temperature of the remainder of the system at the same point on the FFS trajectory. 
Figure [6] shows that the bulk liquid temperature of the total system computed with the 
collected MD snapshots at each interface i is effectively constant and equal to the imposed 
temperature of thermostat. 

Figure El shows that a successful nucleation is likely to have started in a region where 
the local temperature is higher than the temperature of the bulk liquid. In other words, a 
growing bubble is, at least initially, hotter than the surrounding liquid. From interface 6, 
the local temperatures of the largest bubbles is not significantly different from the bulk tem- 
perature. By analyzing the overall volume of the system, we find that it barely changes from 
interfaces ~ 4, whereas it starts increasing from interface 5 (see Figure [T]). Interestingly, 
the system expansion starts at the point where the bubble has cooled to the temperature of 
the surroundings. 

Between interface 5 and the critical (iso-committor) surface 12, there appears to be 
a difference between bubbles that continue to grow and those that shrink. On average, 
the temperature of growing bubbles appears lower than the ambient temperature, whilst 
is not significantly different from ambient. Note, however, the statistical errors are large: 
individual points may not conform to this general trend. It should be noted that the effect 
of local temperature on nucleation may, in practice, be somewhat larger than indicated here 
because the system that we study is always weakly thermostatted. Hence, to return to 
the example of the collapsing bubble, the final temperature reached on collapse would be 
somewhat less in our simulations than in reality. 

The main conclusion of the present section is that local temperature fluctuation play 
an important role in bubble nucleation. A description such as CNT that is based on the 
assumption of local thermal equilibrium cannot capture this phenomenon adequately. In 
the wider context of activated processes in condensed media it is interesting to find an 
example of a rare event where the kinetic energy appears to be a relevant quantity to 
determine the propensity for the subsequent barrier crossing. Such an observation would 
not be surprising for barrier crossings that are largely ballistic, but it has, to our knowledge, 
not been observed in the case of more diffusive barrier crossings. We stress, however, that 
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the effect of local temperature fluctuations on bubble nucleation is not surprising: this 
phenomenon is responsible for the operation of the bubble chamber to detect elementary 
particles. 



4. CONCLUSIONS 

Combining Molecular Dynamics simulations with Forward-Flux Sampling, we have stud- 
ied the kinetics of bubble nucleation in a truncated force-shifted Lennard- Jones liquid. In 
our simulations, we used the volume of the largest bubble as the order parameter to follow 
the phase transition. Using FFS, we computed the bubble nucleation rate and compared it 
with the one predicted by CNT. The simulated nucleation rates appear significantly larger 
than those predicted by CNT. However, the CNT estimates are subject to considerable un- 
certainty, as small statistical errors in the calculated value of the surface tension have huge 
effects on the estimates for the nucleation rate. We note that the suggestion that CNT un- 
derestimates bubble nucleation rates also follows from a comparison of predictions based on 
CNT and DFT— . Analysis of the nucleation pathway indicates that nucleation takes place 
via the formation of compact bubbles. These bubbles undergo shape fluctuations but are 
mostly compact and fairly spherical. We do not find evidence for the existence of ramified 
and percolating void structured. We argue that the difference between our findings and 
those of ref. [3^ is not so much due to a difference in the models studied but to the use of a 
global order parameter in ref. l39|. We find that local temperature fluctuations are important 
for the propensity of bubble nucleation. At the beginning of nucleation event, the incipi- 
ent bubble is always significantly hotter than the surrounding bulk liquid. When the total 
volume of the system starts increasing, the bubbles cools down to ambient temperature. 
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Appendix A: Details on the used thermostat and barostat 

In our simulations, we use a thermostat and a barostat to keep the thermodynamic 
conditions of the system constant. However, we wish to work under conditions where the 
temperature and pressure control do not affect the dynamics of the system appreciably on the 
timescale for molecular velocity fluctuations. To select reasonable relaxation times for both 
the thermostat and the barostat, we performed several trial runs at different temperatures 
and pressures before selecting = 0.093 for the thermostat and Tp = 1.390 for the barostat. 
Figure lA-ll illustrates the (kinetic) temperature fluctuations in the fluid about a value of 
(T*) = 0.900 and the (virial) pressure fluctuations in the same fluid about a value of (P*) = 
0.064. 

The mass of the "piston" in the barostat was chosen such that the timescale for volume 
fluctuations corresponded to the time it takes a sound wave to traverse the simulation box 
(approximately 1 ps).— 



Appendix B: TSF-LJ liquid-vapor phase diagram 

Errington et al.— computed the liquid-vapor phase diagram of the TSF-LJ by means of 
MC simulations. In fig. IB- II we show Errington's data for the TSF-Lennard Jones {vc = 2.5a). 
For the sake of comparison, we also reproduce the liquid-vapor phase diagram data of the 
full Lennard Jones calculated by Potoff et al.— and by Smit^^. and the corresponding results 
for the truncated and shifted Lennard- Jones at Vc = 2.5a calculated by Smit^^ and by 
Trokhymchuk et al.^. In the same figure, we also plot our own simulation data for the 
TSF-LJ potential. We find excellent agreement between our own simulation results and the 
data of ref. |49. 

As was already pointed out by Smit®^, truncating and shifting the potential has a large 
effect on the location of the critical temperature. Figure iB^ illustrates this observation. 



Appendix C: liquid- vapor surface tension computed with MD simulations 

To compute the liquid-vapor surface tension, we start by preparing a half-liquid-half- 
vapor system with 16000 particles. After having equilibrated a liquid by means of an NPT 
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simulation at the coexistence temperature and pressure (T* = 0.785, P* = 0.026 or T* = 
0.855, P* = 0.046) , we prepare a liquid-vapor interface by evaporating half of the system. 
The thickness of the liquid slab and the size of simulation box are chosen sufficiently large 
to make finite-size effects unlikely. Once the system of coexisting liquid and vapor has 
reached thermodynamic equilibrium, we calculate the liquid-vapor surface tension in an NVT 



ensemble (see ref. l65l ). We compute the density profile p{z), and the normal and tangential 
pressure profiles, Pn{z) and Pt{z) along the z direction perpendicular with respect to the 
liquid- vapor dividing interface (see Figure IC^ . and use the Irving-Kirkwood approach to 
compute the pressure tensor—. The normal component of the pressure tensor is expressed 
as 

p„(.) = (pW>t.r-i(i:E?^^TT^x»(f^]''(f^)) . (c-1) 

whereas the tangential component is given by 

^ \ i j>i "'"u l%l ^v/ 

where p{z) denotes the density profile along the z direction, ks is the Boltzmann constant, 
T the absolute temperature, A = x Ly the area of the simulation box in the x-y plane, 
the distance between two particles i and j, and U the TSF-LJ internal energy (with cut-off 
r* = 2.5a). The angular brackets denote a canonical average. For computational purposes, 
we slice the simulation box into a large number of slabs in the z direction, where the width 
of each slab is set to Az ~ 0.1 a (corresponding to 1141 slabs at T* = 0.785 and 1160 slabs 
at T* = 0.855 respectively). The surface tension is given by 

1 = 1 r (PN{z)-PT{z))dz . (C-3) 



2 Jo 

The surface tensions that we obtain from MD simulations are 7* = 0.204 ± 0.007 at T* = 
0.785, P* = 0.026 and 7* = 0.098 ± 0.008 at T* = 0.855, P* = 0.046. 
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TABLES 



TABLE 1: Statistics collected from different sets of FFS "bubble-nucleation" trajectories at 
T* = 0.855 and P* = 0.026. {Wh)i denotes the order parameter at each interface i as defined in 
section 12.21 Nf is the total number of MD configurations (including momenta) collected at each 
interface i, irrespective of whether the bubble will continue to grow or it will collapse after interface 
i. is the average number of particles in the vapor phase inside or near the largest bubble. In 
this set, we include all particles that are inside or within 2.5cj from the large bubble, identified 
using the approach described in section [221 From interface 14, bubble nearly always grows without 
collapse. 





i (W,), Nf Nf 


i (W,), Nf Nf 


i (Wt), Nf Nf 


i {Wh)^ Nf Nf 


25 4012 50 


3 92 1795 79 


6 240 725 121 


9 550 334 188 


12 950 79 264 


1 37 3619 60 


4 120 548 98 


7 300 285 142 


10 700 329 203 


13 1200 147 271 


2 57 1304 70 


5 180 512 103 


8 400 320 170 


11 820 374 215 


14 1900 174 289 
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FIGURE CAPTIONS 



Figure 1. Probability distribution of the number of neighbors N for both the hquid and the vapor 
at coexistence T* = 0.785 and P* = 0.026. Two particles are neighbors if their distance is smaller 
than rc = 1.6(7. 

Figure 2. Diffusion coefficients as a function of time in the MD simulations with a Lowe- Andersen 
thermostat. The effect of u* on self-diffusion coefficient of the TSF-LJ fluid is shown. 

Figure 3. MD simulation runs from the first to the second interface {u* = 500) in a TSF-LJ 
system at T* = 0.855 and P* = 0.026. 

Figure 4. Typical bubble-nucleation pathway in a super-heated liquid at T* = 0.855 and P* = 
0.026. The top left frame shows the initial stage of the bubble, while the bottom right snapshot is 
the critical bubble. 

Figure 5. Fluctuations of the order parameter as a function of MD time steps for a generic FFS- 
interface at T* = 0.855 and P* = 0.026. The jumps in Wb are due to the merge and break-up of 
smaller bubbles. 

Figure 6. Temperature at each interface from MD-FFS at T* = 0.855 and P* = 0.026. The green 
line denotes the temperature of the thermostat used in MD. The entire liquid bulk temperature 
is denoted with black solid square. The local temperature of the "both-backward-and-forward" 
largest bubble (blue open triangle) and the "only-forward" largest bubble (red solid circle) are also 
shown. 

Figure 7. System volume at each interface from MD-FFS at T* = 0.855 and P* = 0.026. An 

expansion of the overall volume of the system is observable from interface 5. 
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Figure A-1. Comparison of temperature (left panel) and pressure (rigiit panel) distributions 
between MD simulations and Gaussian distributions peaked at the same average value (continuous 
curves) . 

Figure B-1. T — p liquid- vapor phase diagram. The filled squares denote the present simulation 
results of the TSF-LJ fiuid (rc = 2.5a). Errington's results for the same system are shown as a 
drawn curve. All other (open) symbols refer to simulation data for other variants of the Lennard- 
Jones model (see text). 

Figure C-1. Profiles of p*{z) (top panel), Pj^{z) and Pj'{z) (bottom panel, black-solid and 
red-dashed lines, respectively) of the liquid-vapor coexisting system at T* = 0.785 from MD sim- 
ulations. We verified that the number densities of the liquid and vapor are the same as those of 
a saturated liquid and vapor at imposed temperature (T* = 0.785), and that the normal pressure 
equals the saturated vapor pressure P* = 0.026 at this temperature. 
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